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FKPP FRONTS IN CELLULAR FLOWS: 
THE LARGE-PECLET REGIME 

ALEXANDRA TZELLA* AND JACQUES VANNESTEt 


Abstract. We investigate the propagation of chemical fronts arising in Fisher-Kolmogorov- 
Petrovskii-Piskunov (FKPP) type models in the presence of a steady cellular flow. In the long-time 
limit, a steadily propagating pulsating front is established. Its speed, on which we focus, can be 
obtained by solving an eigenvalue problem closely related to large-deviation theory. We employ 
asymptotic methods to solve this eigenvalue problem in the limit of small molecular diffusivity (large 
Peclet number, Pe ^ 1) and arbitrary reaction rate (arbitrary Damkohler number Da). 

We identify three regimes corresponding to the distinguished limits Da = 0(Pe~^), Da = 
0((logPe)“^) and Da = 0(Pe) and, in each regime, obtain the front speed in terms of a differ¬ 
ent non-trivial function of the relevant combination of Pe and Da. Closed-form expressions for the 
speed, characterised by power-law and logarithmic dependences on Da and Pe and valid in intermedi¬ 
ate regimes, are deduced as limiting cases. Taken together, our asymptotic results provide a complete 
description of the complex dependence of the front speed on Da for Pe 1. They are confirmed 
by numerical solutions of the eigenvalue problem determining the front speed, and illustrated by a 
number of numerical simulations of the advection-diffusion-reaction equation. 

Key words, front propagation, large deviations, cellular flows, homogenization, Hamilton- 
Jacobi, boundary layer, WKB 
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1. Introduction. In a wide variety of environmental and engineering applica¬ 
tions, chemical or biological reactions in fluids propagate in the form of localized, 
strongly inhomogeneous structures associated with reactive fronts [3^ . These are 
usually established as a result of the interaction between molecular diffusion, local 
growth and saturation, but their propagation can be greatly facilitated by advection 
by a flow. There has been a growing interest in analysing this impact of advection on 
the propagation of reactive fronts, as indicated by the large number of experimental, 
and theoretical studies (e.g., [3S1IM1I11IS]) and [481 [71 149) , respectively). 

Much of this work focusses on the effect of incompressible two-dimensional peri¬ 
odic flows, and in particular on the cellular vortex flow. Introduced by this is a 
steady flow with streamfunction 

(1.1) -ip{x, y) = -U sin{x/e) sin{y/e), 

where U is the maximum flow speed and 2'!t£ is the period in both x and y. When the 
system is confined between parallel, impermeable walls at y = Q and iri, as considered 
in this paper, the flow consists of a one-dimensional infinite array of vortices rotating 
in alternating directions (see Fig. O- These vortices are conhned within cells that 
are bounded by a separatrix connecting a network of hyperbolic stagnation points. 

In the absence of advection, the simplest model of front propagation is the FKPP 
model, named after the pioneering work by Fisher m and Kolmogorov, Petrovskii and 
Piskunov [24] . This model describes the evolution of a single constituent that diffuses 
and undergoes a logistic growth, leading to the formation of a steadily travelling 
front. In the presence of a cellular flow (or more general steady periodic flows). 
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Fig. 1.1. Schematic of the streamlines for the cellular vortex flow with streamfunction {OP. 
confined between parallel, impermeable walls. The half-cells with anticlockwise (clockwise) circula¬ 
tion are denoted by the -\- (—) signs. 


the corresponding advection-diffusion-reaction model admits pulsating front solutions 
that change periodically with respect to time as they travel [5]. 

The behaviour of these pulsating fronts depends on two non-dimensional param¬ 
eters: the Damkohler and Peclet numbers, 

Da = £/(C/r) and Fe = U£jK, 

where t is the reaction time and n the molecular diffusivity, which measure the 
strength of advection relative to reaction and to diffusion, respectively. In the inter¬ 
pretation of our results, we will consider a fixed geometry and a fixed flow, in which 
case the values of Pe and Da are controlled by k and r, respectively. In practice, 
however, it is easier to achieve this control by varying £ and U (see e.g. [35]). 

This paper focusses on the limit of large Pe, relevant to many applications where 
advection dominates over diffusion. This is a singular limit, of course, since the weak 
diffusion leads to the creation of spatial scales that are vanishingly small as Pe —)■ oo. 
These small scales are apparent in Figure |1.2| which illustrates the dependence of 
the front structure on the reaction time by showing snapshots of the concentration 
for different Damkohler number Da at fixed (large) Pe = 250. For small Da (slow 
reaction. Fig. |1.2[ a)), the front spreads across several cells, with high concentrations 
within boundary layers surrounding the separatrix. For intermediate Da (Fig. |1.2[ b)), 
the front is narrower: its leading edge is confined around the separatrix as it invades 
successive cells. For large Da (fast reaction, Fig. die)), the front is very sharp with 
a leading edge that penetrates into the cell interiors. 

The main quantitative characteristic of the front is its long-time speed, c. This 
speed is a function of Da and Pe only when the initial conditions are sufficiently close 
to a step function. Assuming this, Freidlin and Gartner m showed that c can be 
deduced from the principal eigenvalue of a certain linear operator. This eigenvalue 
can be interpreted in the framework of large-deviation theory: specifically, it is the 
Legendre dual of the rate function g{c) associated with the probability density function 
for the position of fluid particles that have been displaced - by advection and diffusion 
- to a distance ct in a time t ^ 1. Intuitively, these particles control the concentration 
near the leading edge of the front which, by linearisation, is approximately of the form 
exp{—t{g{x/t) — Da)), whence the front speed c = x/t = (;“^(Da) is obtained. An 
alternative approach, based on the minimum speed of propagation, leads to the same 
eigenvalue problem, as established in [47119]. 

The eigenvalue problem does not provide an explicit analytical expression for the 
front speed but needs to be solved numerically, through computations that become 
increasingly intensive as Pe —> oo or Da —> oo. In the present paper, we carry out a 
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(a) Regime I 




(b) Regime II 





(c) Regime III 
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Fig. 1.2. (Color online). Snapshots of the concentration 9 for Pe = 250, illustrating the 
different type of fronts depending on the reaction rate. The reaction rate is determined by: (a) 
Da = 4 X 10“^ (corresponding to a front speed c Ri 0.15), (b) Da = 4 X 10“^ (c 0.44j and (c) 
Da = 4 (c Ri 0.67 ). In each case, three successive snapshots separated by a time interval 1:11 {2c) are 
shown. 


detailed asymptotic analysis of the eigenvalue problem for Pe ^ 1 and arbitrary Da. 
This provides simpler, and in some cases completely explicit, expressions for the front 
speed, extracting the dominant scalings and elucidating the physical mechanisms of 
propagation depending on the relative values of Pe and Da. 

Partial results of this type have been derived for slow reaction i.e. for Da = 
0(Pe“^): the dimensionless front speed was argued to scale like c/U = in 

. This scaling prediction is in agreement with rigorous bounds obtained in |31j and 
was confirmed by numerical simulations mms]. It is consistent with the closed- 
form prediction obtained using a homogenization technique which is however only 







4 


A. TZELLA AND J. VANNESTE 


Table 2.1 

The three distinguished scalings of Da appearing in the asymptotics of the front speed c for 
Pe ^ 1. The scalings are associated with three regimes that correspond to the three types of fronts 
depicted in Figure . 2\ In each regime, the speed of the front is expressed in terms of a non-trivial 
function ^i,i = 1, 2, 3, that involves a distinct combination of Pe and Da. The range of validity of 
each expression is also indicated. 


Regime 

Da 

clU 

Range of validity 

I 

O(Pe-i) 

Pe-3/'‘^l(PeDa) 

Da <IC (logPe)“^ 

II 

0((logPe)-i) 

(log Pe)-l%(DalogPe) 

Pe“^ < Da < Pe 

III 

0(Pe) 

'^3(Da/Pe) 

Da 3> (logPe)“^ 


valid for Da Pe In this regime, c is found to be proportional to the square root 
of the effective diffusivity deduced from a linear cell problem [n m m iMi HO] and 
determined in [111113137]. In the opposite limit of fast reaction, i.e. for Da = 0(Pe), 
c can be deduced from the homogenization of a Hamilton-Jacobi equation iniiiHiiis] 
and computed by minimizing a certain action functional |45j . The present paper 
extends these results to provide a complete description of the asymptotics of c as 
Pe —)• oo. 


2. Main results and outline. We carry out an asymptotic analysis of the eigen¬ 
value problem determining c and identify three distinguished regimes, characterised 
by the value of Da relative to Pe. These three regimes correspond to the three types 
of fronts depicted in Figure [0| In each regime, we obtain the front speed in terms 


of a non-trivial function of a combination of Pe and Da (see Table 2.1). The function 
relevant to each regime is obtained by solving one-dimensional problems numerically. 
We moreover show that the three regimes overlap for intermediate values of Da, thus 
confirming that our results cover the whole range of Da. 

Our derivation of c in the first two regimes exploits the matched-asymptotics 
analysis recently carried out by Haynes and Vanneste m- Their paper considers the 
dispersion of particles in an unbounded cellular flow and derives the rate function g 
from which we infer c (after some adaptation to account for the walls). The analysis 
captures the behaviour of the concentration in the interior of the cells at the leading 
edge of the front. In Regime I, the concentration is found to be nearly constant along 
the streamlines (see Fig. |1.2K a)), while in Regime H the concentration is vanishing 
inside the cell’s interior (see Fig. |1.2K b)). In both regimes, a boundary layer around 
the separatrix is crucial for the front dynamics. In Regime HI, where the reaction is 
fast, we rely on a Wentzel-Kramers-Brillouin-Jeffreys (WKB) approach which shows 
that c is controlled by a single action-minimising trajectory [45) . 

We note that our predictions are formal, involving no rigorous estimates of the 
associated errors. Instead, they are verified against values of c derived from (i) nu¬ 
merical solutions of the principal eigenvalue problem, and (ii) direct numerical sim¬ 
ulations of the FKPP advection-diffusion-reaction equation. Figure |2.1| shows that 
the asymptotic expressions for c are in excellent agreement with the corresponding 
values obtained from the eigenvalue problem. 

There are four subregimes in which the asymptotic expressions for c reduce to 
closed forms. The reduced expressions, which in fact cover most of the (Da, Pe)-plane 
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(a) (b) 

Fig. 2.1. (Color online), (a) Comparison between numerical and asymptotic predictions for the 
front speed c as a function of Da and for different values of Pe. The numerical results (solid, black 
lines) are obtained from f3. 1 0]) by solving the eigenvalue problem numerically. The asymptotic 

results (colored, dashed lines) correspond to three distinguished regimes describing the three types of 
fronts shown in Figure m with the associated predictions reported in Table {2^ (h) Same as (a) 
but focussing on small values of Da. 


for Pe ^ 1, are summarised in Table [2^ They are used to verify the overlap between 
regimes mentioned above. In these subregimes, c behaves qualitatively as follows. For 
Da ^ Pe~^, the diffusive approximation obtained from classical homogenisation the¬ 
ory is recovered. For Pe“^ ^ Da <C (logPe)“^, c is proportional to Da^'^'*(logPe)”^/"^ 
and is controlled by the dynamics along the separatrix, with the hyperbolic stagnation 
points at the cell corners playing a negligible role. The range (logPe)“^ <C Da <C Pe 
captures the slow growth of c with Da which, in contrast, can be attributed to the stag¬ 
nation points. The expression for c in this range can in fact be crudely approximated 
as the Da-independent c ~ 7 r/logPe. This is qualitatively similar to the expression 
obtained in miD] using a heuristic approach based on an alternative model, the so- 
called G-equation (see also m for a more rigorous analysis). To our knowledge, no 
equivalent expression has previously been derived from the eigenvalue problem. For 
Da 3> Pe, the reaction is so strong that advection contributes only a small correction 
to the well-known FKPP speed Cq = 217-\/Da/Pe = 

The paper is structured as follows. In section we give a brief derivation of the 
eigenvalue problem for the front speed c. The relation between the eigenvalue problem 
and large-deviation theory is also described there. Sections EE and are devoted 
to each of the three distinguished asymptotic regimes. The explicit expressions for 
c in the four subregimes reported in Table |2.2| are also derived in these sections. 
Comparisons with numerical results are presented in section The paper ends with 
the concluding sectionj^ Technical details are relegated to three Appendices. A word 
of caution about our notation may be necessary: to avoid a proliferation of symbols, 
we use the same letters to denote quantities that are scaled differently in each of the 
three Regimes. Specifically, we systematically denote by /o and q the (leading-order) 
Legendre duals to g and c suitably scaled in each Regime, and by 7 the combination 
of Da and Pe on which c depends transcendentally (this is the argument of each of 
the functions in Table 2.1). This should not lead to confusion since these scaled 
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Table 2.2 

The four subregimes and the corresponding closed-form expressions for the front speed. Here 
V Ri 0.53 ,and Wp denotes the principal real branch of the Lambert W function US- 


Subregime 

Range of validity 

cfU 

la 

Da <IC Pe“^ 


Ib/IIa 

Pe“^ <IC Da <IC (logPe)“^ 

7ry^O(4/3)3/4p^3/4Qjjg Pe)-F4 

Ilb/IIIa 

(log Pe) “ ^ <C Da <C Pe 

7r/Wp(8Pe/Da) 

Illb 

Da 3> Pe 

2'\/Da/Pe(l + 3Pe/(16Da)) 


quantities are used exclusively and independently in each of the sections |4]^ 

3. Eigenvalue problem for the front speed. We investigate the propaga¬ 
tion of a reactive front that is established in the cellular flow with streamfunction 
O. The governing equation is the FKPP advection-diffusion-reaction equation 
that describes the evolution of the reactive concentration 9{x,t). Taking ^ as refer¬ 
ence length and the advective time scale ijU as reference time, this equation takes 
the non-dimensional form 

(3.1) dtO+ u-Ve = + D&r{B), 

where u = {ui,U 2 ) = {—dy'ip,dx'ip) and 

(3.2) ijj{x,y) = — sinxsmy, 

are the dimensionless velocity and streamfunction. Here, the reaction term is r(6) = 
9{1 — 9) or, more generally, any function r{9) that satisfies r(0) = r(l) = 0 with 
r{9) > 0 for 9 S (0,1), r{9) < 0 for 0 ^ [0,1] and r'(0) = supg^g^^ r(0)/0 = 1. 
We take the domain to be an infinite two-dimensional strip with no-flux boundary 
conditions 


(3.3) 


dy9 = 0 


at ?/ = 0, TT, 


and 0 —>■ 1 as X —>■ —oo, 0 —>■ 0 as x —>■ oo, so that the front advances rightwards. As 
initial condition we take 0(x,y, 0) = 0(—x), where 0 is the Heaviside step function. 
Note that our non-dimensionalisation implies that the front speed c will from now on 


be expressed relative to the flow velocity U, as reported in Tables 2.1 and 2.2 


Gartner and Friedlin m showed that the long-time speed of propagation of the 
front can be determined by the behaviour of the solution near the front’s leading edge. 
There, 0^1 and r(0) « r'(O)0 = 0 so that equation (3.1) becomes 


(3.4) 


5*0 


V0 = Pe"M0-|-Da0. 


For t ^ 1, the solution can be written as the multiscale expansion 

(3.5) 9{x,t) = t-l/2gdDa-g(C)) ^i{x, H-) , 

where 


(3.6) 


i = xlt = 0 ( 1 ) 
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is treated as a slow parameter. The Pe-dependent function g{^) is independent of Da 
and characterises the dispersion of purely passive particles. It can be recognised as 
the rate (or Cramer) function of large-deviation theory, which quantifies the rough 
asymptotics of the probability density function of the particle positions for t ^ 1 
|20| . The functions (pi, i = 0,1,2, ■■ ■ , are periodic in x: pi(x -I- 27r, y) = (pi{x, y). The 


boundary conditions (3.3) further imply that 


(3.7) 


dyCpi = 0 at 2 / = 0, TT. 


Substituting (3.5) into (3.4) and equating powers of t ^ yields, at leading order 


an eigenvalue problem for cpQ. Dropping the subscript 0 for convenience, this reads 

(3.8) Pe“^ — u ■ V(p — qdxp + {uiq + Pe~^q^) (p = f{q)(p, 

where q = g'{^) can be treated as a parameter and f{q) = ^g'{0 ~ diO is the 
eigenvalue. The relevant eigenvalue is the principal eigenvalue (that with maximum 
real part) because it corresponds to the slowest decaying solution of ( |3.5| ). The Krein- 
Rutman theorem implies that this eigenvalue is unique, real and isolated, with a 
positive associated eigenfunction ^ > 0. Moreover, /(g) > 0 and is convex [5], so that 
/(g) and g(c) are related by a Legendre transform 


(3.9) 


g(/)=sup(g/-/(g)) 


and /(g) =sup(g/-g(/)). 


With g(/) determined, the front speed may be obtained heuristically by observing 


that the solution to (3.4) must neither grow nor decay exponentially with time in a 


reference frame moving with the front, i.e., for / = xjt = c. This happens precisely 
when g{c) = Da which suggests that the front speed satisfies 


c = g (Da). 


(3.10) 

(Note that subdominant terms in expansion ( |3.5[ ) do not influence the above expression 
for the long-time speed value.) The rigorous treatment in m confirms this to be the 
correct speed. An alternative argument seeks solution to (3.4) of the form exp(—gx-l- 
(/(?) + Da)t)^, recovering the eigenvalue problem (3.8). The front speed is then 
determined from the minimum speed condition 


(3.11) 


. f /( 9 ) +Da 
c = inf -, 

q>0 q 


first introduced in m and easily checked to be equivalent t o (|3.10| ) (see also Ch. 7 in 
[m, m and [471 m). In what follows, we rely on the form ( |3.10| ) of the front speed: 
this makes direct contact with recent large-deviation results obtained in [201 [2T] for 
the problem of a non-reacting passive scalar (i.e.. Da = 0) in an unbounded cellular 
flow which we use in our treatment of Regimes I and II. 


The eigenvalue problem (3.8) - in fact a family of eigenvalue problems paramerized 


by g - plays a central role in this paper. In the absence of flow, /(g) = g^/Pe, 
recovering the classical formula for the speed cq = 2-\/Da/Pe = 2^KjT. For general 
It yf 0, the eigenvalue problem (3.8) cannot be solved analytically. Numerically, it 


can be obtained by straightforward discretisation. Computations are simplified by 
observing that the principal eigenfunction inherits the alternating symmetry of the 


streamfunction (3.2) to satisfy 


( 3 . 12 ) 


(p{x + Tr,y) = (p{x,TT - y). 
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(a) 


(b) 


Fig. 3.1. (Color online), (a) Rate function g calculated numerically for Pe = 250 using Jj3.8^ 
and Jg.Pp (solid, black line). The diffusive approximation obtained from is also shown (dashed 

line). This approximation is equivalent to that derived using homogenization theory and is clearly 
only valid for c 1. (b) Same as (a) but focussing on small values of c > 0. 


Figure [3J] shows an instance of g{c) (here for Pe = 250) obtained numerically by 
computing f{q) on a grid in q, then Legendre transforming (the numerical method 
is described in section]^. Clearly, (7 is a non-trivial function of c, only well approxi¬ 
mated by a quadratic function - corresponding to a diffusive approximation - in the 
immediate vicinity of c = 0. We derive below large-Pe expressions for g that cover 
the entire range of c and, correspondingly, expressions for the speed c that cover the 
entire range of Da. This requires to analyse three distinguished regimes defined by 
distinct distinguished scalings of q, c and Da. 

4. Regime I: Da = 0(Pe“^). The first regime encompasses the limit of Da —?► 0 
which is usually tackled using homogenization theory (see e.g. [IIITIES]). Homog¬ 
enization approximates the advection-diffusion equation for a passive scalar by a 
diffusion equation, in which an effective diffusivity Neff replaces molecular diffusivity. 
This approximation assumes that x = for t ^ 1 and implies that 


(4.1) 


5(c) 


-PeKggC^ and /(g) - Pe Veff 5 ^ 


for c <C 1 and g <C 1 (see (3.6)). For Pe ^ 1, the effective diffusivity for the cellular 

flow [HI Sni 133111] is 


(4.2) 




2i/Pe^^^, with iz « 0.53, 


and was obtained in closed form in m- Figure |3.1| confirms the validity of this 
approximation and demonstrates its limitation to a very small range of c. 

Regime I applies to a broader range of c. It can be analysed following [2T] by 
introducing the rescaling 

g = Pe“^/'‘ 


(4.3) 


5, 


as suggested by the form (4.1) of /(g) as g 
are then expanded according to 

(4.4a) /(g) = Pe" Vo(5) + 0(Pe-"/"), 

(4.4b) (/) = </o + Pe-^/Vi+Pe”^/"</2 


where g = 0 ( 1 ), 

0. The eigenvalue and eigenfunction 


Pe 


- 3 / 4 , 


(/3 + Pe"V4 + 0(Pe-®/4). 
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It is convenient to use the value of the streamfunction ip and the arclength s along 
streamlines as coordinates alternative to {x,y). Substituting (4.4) into (3.8) and using 
- 1 .. 


that dgX = 


we obtain the sequence of problems 


(4.5a) ds (po = 0, 

(4.5b) ds4>i = qdsX(j)i-i, fc = l,2,3, 

(4.5c) \\u\\~'^A(po - dstpi + qdsx4>3 = ||m||“V o<?^’o- 


It follows that 00 = 00 (V') is constant along streamlines and automatically satisfies 
condition (3.12). The functions (pi for z = 1, 2, 3 are polynomials in x{ip, s) of degree 


i with 0-dependent coefficients. They do not satisfy (3.7) and (3.12), but these are 


restored through boundary layers at a; = 0, tt and j/ = 0, tt which we treat below. 


Integrating (4.5c) around a streamline leads to the solvability condition 


(4.6a) 


d 

d0 


a(0) 


d0o 

d0 


= /o0o &(0)- 


In this equation, derived using that ^ i(^l|V0|| ds = A0|jV0|| ^ds [35], a(0) and 

6(0) are the circulation and period of orbiting motion along the streamline 0; they 
are given explicitly by 

(4.6b) a(0) = 8(E'(0) — 0^K'(0)) and 6(0) = 4K'(0), 


where K' and E' are the complete elliptic integrals of the first and second kind [T]. 
Note that (4.6) is analogous to an effective diffusion equation obtained by averaging 


Equation (4.6) can be integrated from the centres 0 = =Fl of the half-cells out¬ 
wards. Here we need to distinguish two types of half cells: the ‘-I-’ half-cells, rotating 
counterclockwise with 0 = — 1 at their centre and exemplified by (x, y) € [0, tt] x [0, tt]; 
and the ’ half-cells, rotating clockwise with 0 = 1 at their centre and exemplified 
by {x,y) G [tt, 27r] x [0,7r]. Using systematically the upper (lower) signs for ‘-I-’ (‘—’) 
half-cells, we write the boundary conditions at the centre as 


(4.7) 00 = 1 and = at 0 = =Fl- 

The first condition fixes an arbitrary normalisation for 0o (because ( |4.6[ ) is linear); 
the second ensures that 0o remains bounded as 0 —>■ =Fl (see [H] for details). The 
solution for 0 —)• 0 determines the Dirichlet-to-Neuman map =^(/o), defined as 


(4.8) 


lim 

i/;— 



±^(/o). 


Since ^{fo) 0 0, 0o has a discontinuous first derivative across the separatrix 0 = 0. 
This is resolved by a boundary layer which we examine next. 

Inside the boundary layer, we use the rescaled variables introduced by m, 


(4.9) 


0 = =FPe^/^0 and a 



IIV0II ds. 


where C is a rescaled streamfunction whose sign is chosen so that 0 > 0 in the interior 
of the ± half-cells. Note that 0 < tr < 8 and that the cell corners correspond to 
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(T = 0, 2, 4, 6 . We denote by th e eige nfunction in the boundary layer, and 

expand this in powers of Pe“^^^ as in (4.4b). To leading order <i)o is a constant, 


matching the interior solution: $o = The higher-order terms * = 1,2 

satisfy forced heat equations, with a the time-like variable. Solving these in exactly 
the manner used in the computation of Neff [m uni [sum] leads to the boundary-layer 
counterpart of (4.8), namely 


(4.10) 


r our ^2 

hm —— = 0 and hm <!>„ =-^9 , 

C^oo dC c-foo ^ d( 4 


A derivation is sketched in Appendix]^ The matching of the derivative of p is ensured 
to leading order provided that 


(4.11) 


hm (pQ —— = =F hm $0 . 

^_>o± “ d-ip c-foo dC 


Equating the right-hand sides of (4.8) and (4.10) then yields 
(4.12) 


where denotes the inverse of Recalling the scaling /(g) ~ Pe“^/o(g), the 
above expression gives the asymptotic form of /(g) in Regime I. Note that this ex¬ 
pression is the same as that obtained previously in |21j for an unbounded domain: 
the difference in boundary conditions arising from the presence of walls at y = 0 , tt 
turns out to be unimportant in this regime. 


The front speed is now determined using (3.10). From (4.4) and (4.12), we deduce 


that 


(4.13) 


g{c) = Pe"^^i(Pe3/^c) -f 0 (Pe"®/^) 


where is the Legendre transform of ^ Solving (3.10) then gives 


(4.14) 


Pe"^/Vi( 7 ) for 7 = DaPe = 0(l), 


where Note that, although this expression is derived assuming formally 

that 7 = 0(1), it will become clear from our analysis of Regime II below that it 
applies for the larger range 7 ^ Pe(logPe)“^. 


Eq. (4.14) shows that for a fixed value of 7 , and thus for constant front thickness 
(since in the absence of advection, this thickness is (kt)^/^ = c oc Pe“^^^, 

which explains the power law that was previously conjectured in [5] and observed in 
the numerical work of |46| . It is also consistent with the rigorous upper and lower 
bounds scaling as Pe“^/^ obtained in m under the assumption that 7 = 0(1). It is 
straightforward to determine '^1 numerically and thus obtain an approximation for c. 
We first calculate for gridded values of /o using s tand ard second-order finite 

differences to discretize (4.6) with boundary conditions (4.7). Inverting gives 


then, by Legendre transforming, Another inversion finally yields '^ 1 . The result 
is shown in Figure |4.1| This demonstrates that is a non-trivial function of its 
argument, implying that a power-law approximation is only valid locally. 
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Fig. 4.1. Large-Pe prediction 4-^4 ^ the speed c scaled by Pe^^^ as a function 0/7 = PeDa. 
The function ^ 1 ( 7 ) in f4-14\ l derived by inverting the Legendre transform of fo (solid black line) 
which is obtained from by solving an ODE ( (4-^ with boundary conditions \4- )■ small-'j 

approximation . I g[ ) (dot-dashed line) is also shown, along with the large-'y approximation obtained 
(i) by solving f4DS\ 

(upper dashed line) 


numerically (lower dashed line) and (ii) using the cruder approximation (4 ■ 20\ 


Asymptotic limits. We now derive two approximations for that result in 
two asymptotic subregimes la and Ib of Table [2?^ Both approximations are based on 
the asymptotic form of ^(/o) that [5T] derived for small and large values of fo ^ Pe/. 
The first approximation uses that 


(4.15) 


'^(/o) —-^/o + 0 (/o) as fo 


0 . 


Introducing into (4.12 1 recovers the quadratic approximation (4.1) for fo{q)- We 


employ ( |4.13[ ) and ( |4.14[ ) to deduce that 

( 81 ^ 7 )^/^ as 7 —^ 0 . 


(4.16) 


■^1(7) 


Figure 4.1 confirms the validity of this approximation. Eq. (4.14) then gives the front 
speed as 

(4.17) 


c ~ ( 8 i/)^/^Da^/^Pe"^/'‘ for Da < Pe 


-1 


The validity of this approximation was previously established in [SHiisn] where it was 
shown that in the limit of Da —^ 0, the front speed is calculated from the quadratic 


approximation (4.1). 


The second approximation uses that 


(4.18) 


^(/o) = 


V2X 


1 + 


log A 


+ 0 ((logA) ) as fo^oo, 


where A is the solution of A^ = 4/ologA and /i « 0.81. Figure 4.1 shows that the 


corresponding approximation for - obtained by numerical evaluation of (4.18), in¬ 


version and Legendre transform - is very accurate when its argument is sufficiently 
large. We emphasise that this approximation, although it requires numerical compu¬ 


tations, is much simpler than (4.12) in that it requires only the solution of algebraic 


equations instead of the solution of a differential equation. A closed-form expression 
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is deduced by solving the transcendental equation defining A asymptotically to obtain 
the leading-order approximation 


(4.19) 


^(/o) 


^(/olog/o)^/^ as fo 


oo, 


noting that the second term in (4.18) is subdominant. This approximation is crude 
because it ignores terms that are O((log/o)“^) relative to the term retained. It 
is n oneth eless useful because it leads to an explicit expression for the speed: us¬ 
ing (4.12) gives that /olog/o and hence, to leading order, that /o ^ 

TT'^v'^q'^/{IQlog q) as g —>■ oo. Ultimately, using f{q) ^ Pe“^/o(g), this gives 


(4.20) 


‘^1(7) ~ (4/3)^^^ 7^/'^(log7) as 7 


00. 


This expression captures the asymptotic behaviour of ‘^1(7) but, as Figure 4.1 shows 


the logarithmic corrections that it neglects are substantially large for finite 7. Using 


(4.14), we deduce the approximation 


(4.21) 


(4/3)^^“* Da^^‘*(logPe) for Pe ^ ^ Da (log Pe) 


where the upper bound corresponds to the distinguished limit of Da associated with 
Regime II. Note that we have dropped a term in log Da using that Pe Da and 
Pe Da“^. Expression (4.21) will be used below to verify the matching between 


regimes I and II. 

5. Regime II: Da = 0(l/logPe). This second regime applies to values of 
Da larger than in Regime I which it continues smoothly. The analysis, which again 
involves boundary layers, is similar to that carried out in |2Ij for the non-reacting 
problem. There are however major differences stemming from the bounded domain 
that we consider; we therefore describe the analysis in some detail. 

Motivated by the observation that / = 0{q^) when q ^ Pe”^'^^ (up to logarithmic 
terms, see the discussion preceding (4.20)), we assume that q = 0(1) and expand the 
eigenvalue and eigenfunction as 


(5.1) 


/(<?)=/o(g) + 0(Pe-i/^) and 


Introducing (|5.1[) into the eigenvalue equation (3.8), we find that the interior solution 


vanishes at leading order: (j)Q = 0. Thus the solution is entirely determined by the 
behaviour in the boundary layer around the separatrix, as the numerical simulations 
hint (see Figure Bb)). 

The boundary layer has a thickness 0(Pe ^^^), as in Regime I; inside, the leading- 
order solution satisfies 


(5.2) 


die 


$0 - 


U 


where $0 is expressed in terms of the rescaled variables (4.9). This can be turned 


into a heat equation along each segment of the boundary layer using the piecewise 
transformation 


(5.3a) 


$ = exp (—gx —/o Fl(cr)) $0, where H{a) =— 


1-2 


'2k/2j 


dcr^ 
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which reduces (5.2) to 


(5.3b) 




This transformation breaks down near the cell corners where ||it|j vanishes. There, 


different rescaled variables, namely {X,Y) = Pe^^‘^{x,y), are required to solve (3.8). 
The solution that is obtained to leading order, n amel y $o = X~^°^{XY) for some 
function <i>, can be matched with the solution of (5.3) upstream and downstream of 


the corner. This leads to jump conditions at each corner reading 


(5.4) lim $((, cr) = (16Pe) lim $(^, cr), for fc = 0, 2, 4, 6, 

CT—cr— 


(see [H] for details). Combining these jump conditions with (i) the relation between 
<I> downstream of each corner and upstream of the next corner that follows from 
(5.3b), and (ii) the symmetry (3.12) and boundary conditions (3.7) results in the 
eigenvalue problem 


(5.5) 


(16Pe)*/2$(C) = (/C$)(C) 


(see Appendix]^. Here $ is a vector grouping the four solutions downstream of each 
corner, that is, $(C,cr) for a = O'*", 2+, 4+, 6+, and /C is a 4 x 4 matrix operator 
that depends explicitly on q and /q. Its entries are linear combinations of the linear 
integral operators ‘K^ defined by 

(5.6) (j£±cj,)(^) = / e-(‘;^«')'/8<i>(C')dC', 

VStt Jo 


for an arbitrary function $. 

An expression for /o is now obtained by considering the principal eigenvalue of 
K. Let A denote this eigenvalue. Introducing into (5.5) and solving for /o gives 


(5.7) 


/o = 


2 log A 
log(16Pe) ’ 


where A = A(g,/o). 


Note that even though log 16 provides an asymptotically negligible correction to 
log Pe, it turns out to be significant for the large-but-finite values of Pe we consider 
and is therefore better retained. 

Equation (|5.7[) is transcendental. It is solved numerically by first discretising K. 


to find A(g, /o) as the eigenvalue of a matrix, then solving (5.7) iteratively, using the 
straightforward scheme 


(5.8) 


f(^) _ 


! log A (g, 
log(16Pe) 


n = 1, 2, 


taking = 0 as initial guess. This guess is reasonable when g <C 1 in which case 
/o <C 1. As the value of g increases, the sequence of corrections generated by (5.8) 
become increasingly important, and increasingly larger values of Pe are needed for 
the leading-order approximation to be accurate. 


The front speed can be derived from the solution /o = /o(g,Pe) to (5.7) by 
Legendre transforming with respect to g to obtain g(c), then solving g(c) = Da. 
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This leads to c as a transcendental function of Da and log(16Pe) that can approx¬ 
imated numerically, starting with the estimate for /g obtained by iterating (5.8). 


This approach does not make explicit the scaling relation that characterises Regime 
II, however. To obtain this, we approximate X{q, fo) in (5.7) by A(( 7 , 0), leading to 
/o ~ /g^^ = 21ogA((7,0)/log(16Pe), and hence to 


(5.9) 


5(c) 


i^2(log(16Pe)c) 

log(16Pe) 


where denotes the Legendre transforms of 2 log X{q, 0) with respect to q. The front 
speed asymptotics 


(5.10) 


*^2(7) 

log(16Pe) 


for 7 = log(16Pe) Da, 


where ^2 = '^2 follows. We emphasise that this approximation is asymptotically 
consistent for q — 0(1) since /g —>■ 0 as Pe — 00 . As we show shortly, its accuracy 


is poor for finite Pe and the complete solution to (5.7), which treats l/log(16Pe) as 
0 ( 1 ), is preferable. 


Figure |5.1| shows the behaviour of ^2 obtained numerically for a range of values 
of 7 = log(16Pe) Da. The range is limited because the matrix associated with the 
discretised version of K, (with /g = 0 ) becomes ill conditioned when 7 > 1 , leading 
to numerical inaccuracies in the principal eigenvalue A(( 7 ,0). The complete solution 
to (5.7) leads to a (Pe-dependent) approximation to clog(16Pe) which, in contrast, is 


well conditioned over a broad range of 7 ; this approximation is shown in Figure [5T] for 
four values of Pe. The results indicate that the logarithmic corrections included in the 


complete solution are negligible for 7 1 , with (5.10) providing a good approximation 


but significant for larger 7 when they are seen to decrease very slowly as Pe increases. 
The results are also consistent with the behaviour clog(16Pe) ~ ‘^ 2 (c) ~ tt for 7 ;g> 1 
derived below. 

Asymptotic limits. There are two asymptotic approximations of the front speed 
in Regime II, corresponding to 7 ^ 1 and 7 1 and identified as subregimes Ila and 


Ilb in Table 2.2 For the first, we approximate '^2 in (5.10) based on the asymptotic 


form of X{q, 0) for q 1 derived in [2T]. For such q, the jumps in (5.4) are negligible 


and the boundary-layer solution can be expanded in powers of q, whence it is found 
that X{q,0) ~ exp(2/r^) where /i = 7 r^i/g^/ 4 . It follows that 


(5.11) 


'^2(7) 


d/2 


(4/3) 


3/4 ^ 3/4 


for 7 < 1 , 


and, using (5.10), that the front speed is that reported in Table 2.2 This Regime 


Ila asymptotic expression coincides with that found in Regime Ib as (4.21), thus 
confirming the matching between Regimes I and II. 

The second asymptotic approximation corresponds to 7 ^ 1, hence g 3> 1. In this 
limit, the eigenvalue X{q,fo) of K can be derived from a scalar eigenvalue problem 
which we derive and solve asymptotically in Appendix From this solution, we 
deduce the asymptotics (C.9) for /g. Taking the Legendre transform gives g{c) ^ 
8 Pe ce“’^/'^/ 7 r, which we invert to obtain the front speed in Regime lib as 


(5.12) 


Wp( 8 PeDa"^) 


for 


log Pe 


< Da < Pe, 


















^2(7) 


FKPP FRONTS IN CELLULAR FLOWS 


15 




(a) 


(b) 


Fig. 5.1. (a) Large-Pe prediction \5. 1 for the front speed c scaled by log(16Pe) as a function 

of-y = \og(16Pe)Da = 0(1) (lower dashed line). The function ^2 1 0]l is obtained by computing 

the principal eigenvalue of J5.5p for fo = 0 numerically. The small-y approximation (5.11 ) (upper 
dashed line) and large-y approximation *^2 ~ (dashed-dotted line) are also shown. The four thin 
solid black lines correspond to higher-order corrections to (5.10^ obtained for Pe = 50, 125, 250 
and 500, with the arrow pointing in the direction of increasing Pe. (b) Higher-order corrections 
compared to the large-y approximation Il5.12\l (dashed lines). 


where Wp denotes the principal real branch of the Lambert W-function ( solut ion 


of W{z)e^^^'> = z, see [T]) The upper bound of the range of validity of (5.12) is 
determined by comparison with the results in Regime III in the next section. 


Figure 5.1 shows that the approximation (5.12) is very good for Pe = 50, Pe = 125 


and excellent for Pe = 250 and 500. We note that, since Da ^ Pe a nd Pe S> Da , it 
is consistent to approximate Wp( 8 PeDa“^) by logPe [T] to reduce (5.12) to 


(5.13) 


log Pe 


This approximation is poor for finite Pe because of the neglect of logarithmic error 
terms. It is useful in that it shows that both the small-y and large -7 approximations 


lead to the same scaling (5.10) for the front speed, with ‘^ 2 ( 7 ) —tt as 7 —?> 00 . 


We note that an expression qualitatively similar to (5.13) was obtained in mm 
using the so called G-equation, a model alternative to (but not derived from) the 
FKPP model when applied to fast reaction. This expression suggests that the front 
speed c is independent of Da for a range of Da; as the more complete approxima¬ 


tion (5.12) shows and Figure 5.1 confirms, there is in fact a slow growth of c with 


Da. This growth is actually logarithmic, as can be made explicit by improving the 


approximation of (5.12) to include the first-order correction to (5.13) and obtain 


(5.14) 


TT TT log Da 


log Pe log Pe 


6. Regime III: Da = C)(Pe). This final regime corresponds to a fast reaction 
and may be referred to as a geometric-optics regime. Our analysis of Regime lib (and 
specifically, (C.9)) suggests that Regime III emerges for q = 0(Pe). We therefore 
introduce the rescaling 


( 6 . 1 ) 


q = Peq, where g = 0(1), 
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into the eigenvalue problem ( |3.8[ ). We then expand the eigenvalue according to 
(6.2a) fiq) = Pefoiq) + 0(1) 

and assume that the eigenfunction takes the WKB form 


(6.2b) 


= e 


— Peto 


(« + 0(Pe-i)), 


where w and a satisfy the same boundary conditions as (j). Substituting (6.2) into 
(|3.8|) leads to 


(6.3) 


i/(Vw, x) = /o, where H = || Vr(;||^ + u ■ Vw + 2qdxW + uiq + if 


can be regarded as a Hamiltonian. This nonlinear eigenvalue problem has been ob¬ 
tained for general flows by Freidlin, Evans and Souganidis, and Majda and Souganidis 
(see m, [Mj and [28]). It can be interpreted as the cell problem arising in the ho¬ 
mogenisation of the Hamilton-Jacobi equation dtw + ||Vu>|p -I- u ■ Ww = 0 and has 
been shown to have a unique solution /o for each value of q [55|. 

Eq. (6.3) cannot be solved analytically in general, and direct numerical solutions 
are rather involved (see e.g. [23] for the specihc case of the cellular flow). Here we 
exploit a variational formulation which expresses /o, or rather its Legendre dual, the 
rate function go (such that g{c) = Pe50(c) -I- 0(1)), in terms of a minimum-action 
principle. We derive this variational formulation by considering a time-dependent 
version of (6.3), namely the Hamilton-Jacobi equation 

(6.4) 


dtW + H{'Vw,x) = 0, 


noting that we can expect 
(6.5) 


foiq) = - lim 


w(x, t) 


for a wide range of initial conditions w(x^Q). The solution of (6.4) can be written in 
terms of action-minimising paths <p(-) = € K x [0,7r], specifically as 


(6.6) w(x,t) = inf I / L^(9?(s),<p(s))ds 

¥>(•) (Jo 


ip{t) = x> , for t > 0, 


assuming that w{x,0) = 0 (e.g., [IS]). Tbe Lagrangian Lyj is derived by taking the 
Legendre transform of H to find 


(6.7) 
where 

( 6 . 8 ) 


Lw{^{s),^{s)) = L{ip{s), ip(s)) - #i(s). 


L(<p(s),tp(s)) = -||93(s) -ti((p(s))f. 


Using (6.6) and (6.7), we rewrite (6.5) as 

(6.9) /o(g) = - lim ^ inf |g(^i(0) - gx-f / L(<p(s), <p(s)) ds 

t-yoc t ¥>(.) ( Jo 


(p(t) = X 


Without loss of generality we choose x = = 0 and leave y = 4>2 {t) undetermined. 

This is possible because changes to their values lead to 0(1) changes to the infimum 
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and therefore leave /o unaffected. We further make the transformation <p(s) i—>■ —cp{t— 
s). This leaves the Lagrangian ( 6 . 8 ) unchanged and enables us to rewrite /o as 


(6.10) fo{q) = \im ^ sup \ q^i{t) - [ L{(p{s),(p{s)) ds 

t cp{-) I Jo 

We now introduce c = (pi (t) jt to obtain that 

/o(g) = sup (gc- lim ]■ inf | / L((^(s), ¥?(s)) ds 
c V ^ vi-) [Jo 


<^( 0 ) = ( 0 ,-) 


¥?( 0 ) = (0,-),¥’(i) = (ct,-) 


where the dependence on specific values of ^ 2 ( 0 ), ^ 2 {J) is dropped. Recognizing the 
Legendre transform, we obtain the rate function 


( 6 . 11 ) 50 (c) = lim - inf 

t—>-oo t <p{-) 


HJ^is),ip{s))ds 


<^( 0 ) = ( 0 , •),¥’(t) = (ct,-) 


This gives 50 (c) in terms of the action-minimising path - or instanton - ¥’*(•). 

We make three remarks. First, the result (6.11) follows directly from an ap¬ 
plication of the Freidlin-Wentzell (small noise) large-deviation theory (see |18).|171 
Ch. 6 ] and [16]) to the dispersion of passive particles in the flow u. Thus Regime 
III can be regarded as lying at the intersection between large-t large-deviation theory 
as used in this paper, and small-noise (large-Pe) large-deviation theory: that their 
results coincide indicates that the two limits t —)■ 00 and Pe —)■ 00 commute. Second, 
the asymptotics of the principal eigenvalues of a broad class of second-order elliptic 
operators can be obtained using a variational approach [34]; thus, (6.11) could be 
alternatively derived by application of the relevant results in |34| . Third, since ( |6.3| ) 
is the cell problem for the homogenisation of a Hamilton-Jacobi equation [HIH], 

(6.11) provides a variational route to derive the homogenised Hamiltonian /q. 
Computing the right-hand side of (6.11) becomes considerably easier by observing 

that we may take the minimising path to be periodic, in the sense that 

(6.12) (^(s-I-r) = <^(s)-I-(27r, 0), where t = 27r/c. 


Using that L{(p{s), (p{s)) ds = n fj L(cp(s), ¥>(s)) ds reduces (6.11) to 

(6.13) 50 (c) = - inf I f L{ip{s),(p{s))ds 

T V{-) [Jo 


<^(0) = (0,'),¥’(t) = (Stt,-) 


Recalling the scaling 5 (c) ~ Pe 50 (c) and letting cr = cs in the above expression, we 
finally obtain the rate function as 


(6.14) 

where 


5 (c) = Pe^ 3 (c), 


(6.15) ^ 3 (c) = ^ inf I f \\cip'[a) - u{ip{a))fda 

Stt ¥>(•) [Jo 


(p{2Tr) = (p(0) + (27r, 0) 


The front speed in Regime HI follows as 

(6.16) c ^ *^ 3 ( 7 ) for 7 = Da/Pe = 0(1), 
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Fig. 6.1. Large-Pe prediction for the front speed c valid for 7 = DajPe = 0(1). 

^3 is calculated numerically by minimizing (solid black line) and compared with the small-c 

asymptotic approximation (5.12^ (lower dashed line), the large-c asymptotic approximation (6.19)) 
(upper dashed line), and the bare speed cq = 2^/7 (dotted line). The inset focuses on smaller values 
of 1 (after 


where *^3 = . The authors derived this result previously using a different ap¬ 

proach, directly related to the Freidlin-Wentzell small-noise large deviation, that 
bypasses the eigenvalue problem (3.8) [45]. While the present derivation is more 


involved, it highlights the relation with the eigenvalue problem and hence the connec¬ 
tion between the three regimes. 


The minimization problem (6.15) provides an easy way to compute the instanton 


and thus the front speed numerically. Its solution is straightforward to obtain using 
MATLAB’s optimization toolbox. We first start with a large value of c and use a 
standard first-order finite-differences to discretize cr in A^ = 250 equidistant points. 
The resulting discrete action is then minimized using the routine fminunc that is 
seeded with the straight line ¥^*( 5 ) = (cs, 7r/2) as initial guess. We then iterate over 
a range of values of c using the previously determined path as an initial guess to find 
the next minimizer. Figure 2 in [45j shows characteristic examples of instantons tp* (s) 
that are obtained for different values of c. These are close to a straight line when c is 
large and follow closely a streamline near the cell boundaries when c is small. Figure 


6.1 shows the behaviour of c as a function of 7 deduced from (6.16). 


Asymptotic limits. Closed-form expressions for c are derived in [IS] for two 
asymptotic limits, corresponding to 7 <C 1 and 7 1 and referred to as subregimes 

Ilia and Illb in Table [2?^ We sketch the derivation here for completeness. 


For 7 <C 1 and hence c <C 1, the instanton follows a streamline close to the cell 
boundaries, departing from it only for y « 7 r/ 2 . The action (6.15) is minimized when 
= {x{a), y{(T)) satisfies cy' « — cosxsiny (so that the instanton and flow speeds 
differ only in the x-direction). Exploiting symmetry to consider 0 < cr < 7 r /2 only, 
with x(0) = 0, 2 /( 0 ) = x{'k/2) = 7 r /2 and y'{'K/2) = 0, we can divide the instanton path 
into two segments. In region 1, where x <C 1, the integrand in (6.15) is approximately 
(ex' —X cos 2 /)^, leading to the Euler-Lagrange equation c^x" = x (since cy' « — sin 2 /). 
In region 2, ^ 1, ex' = sin x cos 2 / ~ — sinx and cy' = — cos x sin 2 / « —ycosx. 

Matching between the solutions in their common region of validity x, 2 / 1 (the cell 
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corner) gives the approximation 


(6.17) 




{Ci{a),C2{(j)) for cr < 7 r /2 

(^2(^2 - O'), <73(71/2 - cr)) for cr > c 


■where <7i(cr) = 4exp(— 77 /( 20 )) sinh((T/c), < 72 ( 0 -) = 2tan“^(exp(—cr/c))) and C^{(7) = 
4exp(—7r/(2c)) cosh(i7/c). Exp ression (6.17) is in very g ood a greement with our nu¬ 
merical solution. Using (6.17) gives the integrand in (6.15) as {ex' — xcosyY « 
16exp (— tt/c) cosh“^ (o'/C), leading to 


(6.18) 


^ 3 ( 0 ) ~ 4 X {2lT:)ce where c ^ 1 


and the factor 4 appears because, for cr S [0 27 r], the solution (6.17) repeats 4 times, 
up to symmetries. Inverting (6.18) yields 


(6.19) 


IUp(87-i) 


for 7 = Da/Pe 1, 


that is, the same expression as (5.12) found as Regime Ilb. This verifies the matching 
between Regimes II and III. 

The second asymptotic limit corresponds to 7 ^ I, hence c ^ 1. In this case, 
the instanton path is approximately a straight line, with expansion 


( 6 . 20 ) 


c^*(cr) = (cr,j/o)+c ^{xi{a),yi{a)) + 0{c '^) 


where xi, yi are 27r-periodic functions satisfying xi(0) = yi{0) = 0. Substituting 
into (6.15) and minimising with respect to yo, Xi(cr) and yi{a) gives Xi(cr) = 0, 


2/1 (cr) = —2 sin cr sin 7/0 and y^ = 7 r/ 2 , leading to 

(6.21) %(c) = cV4-3/8-bO(c-2). 


Using (6.16) finally leads to the asymptotics of the speed 


( 6 . 22 ) 




for 7 > 1 . 


The leading-order term in (6.22) is the bare speed cq, unsurprisingly since reaction is 


so strong in this regime that advection has a small effect on the front evolution. The 
second term in the expansion is necessary for a good agreement between asymptotic 


and full results (see Fig. 6.1). 


7. Comparison -with numerical results. We compare our predictions for the 
speed c derived in each regime with the corresponding values obtained from (i) the 
numerical evaluation of the principal eigenvalue in ( |3.8| ), and (ii) direct numerical 
simulations of the FKPP equation (3.1) with r{9) = 6{1 — 9). For (i) we use a 


standard second-order finite-difference discretization of (3.8). The resulting matrix 


eigenvalue problem is solved for a range of values of q using MATLAB’s routine eigs. 
We choose the spatial resolution A to satisfy tt/A = 750 in both directions. 

For (ii) we discretize (3.1) using a fractional-step method with a Godunov splitting 
in which we alternate between independent advection, diffusion and reaction steps. 
The advantage of this method is that it is simple and cheap to combine a high- 
resolution finite-volume method for the advection equation dt9 -|- 'U • V 0 = 0 , with an 
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alternating-direction implicit method for the diffusion equation dtO = Pe~^A0, and 
an exact solution of the reaction equation dtO = Dar(0). The advection equation is 
solved using a first-order upwind method that includes a minmod limiter to account 
for second-order corrections (see for more details). This is a stable scheme as 
long as the Courant-Friedrichs-Lewy (CFL) condition is satisfied. We choose the 
spatial resolution A to satisfy tt/A = 400 when Da < 1 and tt/A = 750 otherwise. 
This way we ensure that A/tt > 10“^min((5i,^ 2 ) for all values of Pe and Da where 
(5i = 0(Pe“^^^) and 62 = 0(Pe“^^^Da“^^^) are the characteristic thicknesses of the 
boundary layer and front, respectively. The time-step is controlled by the CFL number 
that we set to be equal to 0.8. 

To make the computational domain finite, we set artificial boundaries at cc = 
ztNir, with N = 15 when Da < 1 and N = 5 otherwise, so that boundary effects are 
negligible. A larger domain is necessary for smaller Da values because the front width 
is larger (see, e.g., Fig. |1.2[ a)). We impose absorbing boundary conditions using a 
zero-order extrapolation at each of the four boundaries. We modify the computational 
domain to track the front for a long time: each time the solution at a; = (N — l)7r 
becomes larger than e = 10“®, we eliminate the nodes with —Ntt ^ x ^ {—N + l)7r 
to the left of the front and add new nodes with Ntt ^ x ^ {N + l)7r to the right of 
the front where we set 9 = 0. The front speed is insensitive to the precise value of e. 
We calculate the speed of the front by considering the left and right endpoints of the 
front, xj{t) and xf(t), defined as 

(7.1) xj{t) = min{x : 9{x,t) = 1 — e} and x'^{t) = max{x : d{x,t) = e}, 


which we determine using a third-order polynomial interpolation. We calculate the 
large-scale speed of the front from a linear fit of xf (t) that we obtain for values of t 
sufficiently large for x+(t) — xj{t) to remain approximately constant. The results are 
not sensitive to the exact value of e: comparison with results obtained for e = 0 . 001 , 
0.01 and 0.1 resulted in less than 1 % of difference in the speed of the front. 

The two sets of numerical results are shown in Figure |7.1| along with the cor¬ 
responding prediction for each regime, respectively derived from (4.14), (5.10) and 
(6.16). The speeds obtained from the eigenvalue equation (3.8) are in excellent agree¬ 
ment with the corres pond ing values obtained from the full numerical simulations of 
the FKPP equation (3.1). This is especially the case when Da > lOPe”^. For 
Da « 5Pe“^, we observe a small dependence of the speed value on the threshold e 
that is used to define the right endpoint of the front (see inset in Figure [7T| ( a)). This 
dependence is due to the particularly long integration times and computational do¬ 
main that are necessary to capture this slowly advancing, wide front (see Fig. |1.2[ a)). 
As Da increases to 0(1) values and beyond, the solutions to (|3.1|) and (3.8) become 


progressively localized, with the smallest lengthscales being 0 ( 62 ), which are challeng¬ 
ing to resolve when Pe ^ 1. This is partly reflected in Figure [7T| (b) where for the high 
values Pe = 250, 500, the agreement between the two sets of numerical results is not as 
close as for the moderate values Pe = 50,125, with the difference increasing with Da. 
In Figure [7T| (c), where the speed is unsealed, the agreement is excellent. However, we 
were not able to obtain sufficiently accurate speed values when Da/Pe = 0(1) from 
either ( |3.1[ ) or ( |3.8[ ) due to the numerical limitations when Da, Pe ^ 1. 

It is clear that in all three regimes, the asymptotic predictions become increasingly 
accurate as the value of Pe increases. In Regime II, the agreement is very good for 
all values of Pe when Dalog(16Pe) is small. However, when Dalog(16Pe) is large, 
we need to employ higher-order corrections to (5.10) (which are obtained via (5.8)). 
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(a) Regime I 



(b) Regime II 



(c) Regime III 



Fig. 7.1. (Color online). Comparison between asymptotic and numerical results of the front 
speed c in all three regimes and for various values of Pe. Solutions of the eigenvalue problem are 
shown in thick, grey (colored) solid lines. Results of the full numerical simulations are shown as 
symbols. T he das hed thin lines are the asymptotic predictions, (a) Regime I show ing the asymptotic 
prediction \4-14 ^ valid for Da = 0{Pe~^); the diffusive approximation {4-11) is also shown in 
the inset that magnifi es the small-PeDa region (dashed-dotted line), (b) Regime II showi ng the 
asymptotic prediction jS. 1 0\l valid for Da = 0((log(16Pe)“^); the corrections to prediction \5. 1 0]) 
obtained by iterating \5.8}) are also shown (thin solid black lines). The inset focuses on two values 
of Pe = 50,500 and on a smaller region of values of Dalog{16Pe). (c) Regime III showing the 
asymptotic prediction valid for Da = 0{Pe), with inset focussing on small values of DafPe. 
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These capture very well the slow growth of the speed values when Dalog(16Pe) ^ 1, 
particularly so for Pe = 250, Pe = 500. In Regime III, the agreement is excellent for 
all values of Pe. 

As expected, the asymptotic expressions are valid over a broad range of values 
of their argument, restricted only by the range of validity of each regime. Taken 
together, they cover the entire range of Da to provide convenient approximations for 
the front speed c, including when Pe and/or Da are so large that direct numerical 
computations are challenging. 

8. Conclusion. In this paper, we study the classic problem of FKPP front prop¬ 
agation in a cellular flow. We examine in detail the asymptotic form of the front 
speed c in the limit of large Peclet number Pe corresponding to a diffusion that is 
weak compared to advection, and for arbitrary values of the Damkohler number Da, 
i.e., arbitrary reaction rate. This is achieved by a careful asymptotic analysis of the 
two-dimensional eigenvalue problem from whose solution c can be deduced. This is 
complicated by the non-uniformity of the problem: depending on the relation between 
Da and Pe, different regimes emerge which require different asymptotic methods and 
lead to different expressions for c. 

Specifically, we identify the three distinguished regimes listed in Table [2T| In each 
regime, the front speed is given in terms of a transcendental function of a suitable 
combination of Pe and Da. Each function is determined by solving a (Pe-independent) 
one-dimensional problem: an ordinary differential equation in Regime I, an integral 
eigenvalue problem in Regime II, and an optimisation problem in Regime III. These 
problems need to be solved numerically in general, though at a much reduced compu¬ 
tational cost compared with the original two-dimensional eigenvalue problem thanks 
to the dimensional reduction, the independence on Pe, and the single scale of the 
solution. Closed-form expressions are obtained by considering asymptotic limits of 
the one-dimensional functions characterising Regimes I, II and III, leading to the 
subregimes listed in Table |2.2| By verifying that the same expressions for c can be 
obtained by suitable limits of both Regimes I and II on the one hand, and of both 
Regimes II and III on the other, we confirm that our formulas cover the full range 
of values of Da. We emphasise that the closed-form formulas valid in the various 
subregimes apply to most of the (Pe, Da)-plane for Pe ^ 1, with the more complex 
distinguished expressions only required in the comparatively narrow regions defined 
by DaPe = 0(1), Da logPe = 0(1) and Da/Pe = 0(1). 

Our analysis reveals previously unchartered behaviour. Only two sublimits are 
intuitively obvious: the first (Illb, Da ^ Pe) arises when the reaction is so fast that 
advection can be neglected, so that the front speed is the familiar bare speed, dimen¬ 
sionally ciiib = Co = 2^njT = 2t7-\/Da/Pe, obtained in the absence of flow. The 
other obvious sublimit (la. Da <C Pe“^) arises when the reaction is slow enough that 
the front spreads across many flow cells; in this case, homogenisation results which 
describe the combined effect of advection and diffusion through an effective diffusivity 
Kefif apply, and the front speed is estimated by replacing k by Neff = 2NPe^^^K in the 
bare speed to obtain cia = . These two explicit expressions provide 

estimates for c for extreme values of Pe, but since their ratio cia/cmb = v^Pe^^^ is 
asymptotically large, they provide little indication (bar a lower bound for cmb) for 
the front speed for Da away from these extremes. Our asymptotic results, in contrast, 
pinpoint the behaviour of c. They describe, in particular, the very slow growth of 
c with Da in Regime Ilb/IIIa where, to the lowest order ignoring logarithmic cor¬ 
rections, c ~ 7r/logPe is independent of Da. This scaling, proposed heuristically in 
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mm, is here derived in two ways, from the integral eigenvalue problem of Regime II 
and from the optimisation approach of Regime III. It can be traced to the behaviour 
of fluid-particle motion near the cell corners: the front in this regime is controlled by 
motion along the separatrix which is fast along most of the separatrix but very slow 
near the corners since these are stagnation points. As a result, the motion of particles 
determining the front is akin to a random walk on the lattice of stagnation points. 
It is not difficult to show that the relevant waiting time, namely the typical time by 
diffusion to move particles across the stagnation point scales like logPe, thus explain¬ 
ing the form of c. A more complex dependence on logPe holds in the entire Regime 
II, reflecting the same physical phenomenon although complicated by a non-trivial 
behaviour between stagnation points. 

We conclude by noting that most of the rigorous work on the asymptotics of 
FKPP front speed focuses on a single large parameter, namely the Peclet number, 
assuming either that Da = 0(Pe“^) <C 1 [HI [Ml [El (SO] or that Da = 0(Pe) 1 

Esns]. Our analysis and numerical work demonstrates the richness of the problem 
when the Damkohler number is allowed instead to take a broad range of value. This 
richness no doubt extends much beyond the specific cellular flow considered in this 
paper; extensions that demonstrate this for a wide class of flows would be desirable. 
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Appendix A. Boundary-layer analysis in Regime I. We establish expres¬ 
sion (4.10). The derivation is essentially identical to the one in m Appendix A.2] 
and is detailed here for completeness. The differences lie in the matching conditions 
(A.2) but, despite these differences, the derivative (4.10) of the eigenfunction remains 
unaltered. 

Inside the boundary layer, we use the rescaled variables and denote solutions 
in the ± half cells by C). The alternating symmetry (3.12) reads 


(A.l) 




This condition, the boundary conditions (3.7) and continuity across the half-cells 
imply that 


(A.2a) 


(A.2b) 


4>±(0,a) = 4>=F(0,a-b4), 
$±(0,a) = $T(0,a), 


dc«’'^(0,cr) = 0 

for 0 < CT < 2 and 4 < cr < 6, 
9^$'^(0,cr) = —d<^$“(0,cr) 
for 2 < cr < 4 and 6 < cr < 8. 


Introducing expansions of the form (4.4) into (3.8) leads to the sequence 
(A.3a) — da) =0 at 0(1), 

(A.3b) _ a,) ^ = 0 at 0(Pe-'=/p, for fc = 1, 2. 


The only admissible solution to (A.3a) is a constant: $q = <I>o = const. Expressing 
||m|| in terms of cr, ( A.3b|) becomes 


(A.4) 


{dcc-9<r) ^k=TF{a)q^t 


1 ’ 
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where F{a) = {2a — a^)~^l'^ for 0 < cr < 2, F{a) = 0 for 2 < cr < 4, and F{a + 4) = 
—F{a). It follows that the fc = 1 solution to (A.4) satisfies cr + 4) = —(^, cr) 

which once combined with (A.l) yields <i)j’^(C,cr) = —cr). The matching condi¬ 
tions ( A.2[ ) now become 

(A.5a) < 1 >)*"( 0 , ct) = —cr), cr) = 0 for 0 < cr < 2 and 4 < cr < 6 , 

(A.5b) cr) = 0, cr) = —cr) for 2 < cr < 4 and 6 < cr < 8 . 


Defining G{a) by G'(cr) = F{a) with G(cr) da = 0, we write the solution to (A.41 
for fc = 1 in terms of G and p as 


(A. 6 ) 


(C, cr) = ±9 (G(o-) + fc(cr, C)) ^’ 0 - 


Here p{C,a) satisfies d^^-p = d^p with p —>• 0 as C —>■ oo. The boundary conditions 
on the cell boundaries impose that d(^p{Q,a) = 0for0<cr<2, 4<cr<6 and 
p(0, cr) = —G(cr) otherwise. The problem describing p is essentially the same as the 
problem solved by m in the Appendix (the exact correspondence is achieved upon 
multiplication of p by — 2 / 7 r and its translation so that cr i—> cr — 2). The key result is 


(A.7) 


nOO 

■/ p(C, 0) dC =-2zc, 
Jo 


where v is the constant defined in (4.2 1 . 

An approximation to and as ^ oo is obtained by noting that the 

leading-order behaviour of 4)^ at large values of C is controlled by its average around 
the streamline so that 


(A. 8 ) 






da 


for fc = 1 , 2 , as C oo. 


We integrate (A.41 first over a, then over ( to obtain that lim^_,,oo = 0 and 

(A.9) 


lim = --q^ 

C—Too 4 


daF(a) 


dCp(Cw) = 


dCp(C,0) = - — 


where we have combined (A.l I and (A.2 1 to find that = 0 for C = Oj fc = 1, 2 

when 0 < a < 2. This derivation uses that p(C, a) d( = —d(p(0,a) = 0 for 

0 < a < 2, that F(a}da = tt, and (A.7). 

Appendix B. Eigenvalue problem in R egime II. We derive the explicit 
form of the asymptotic eigenvalue pro blem (5.5). The functions 4>^, rel ated to the 
eigenfunction 4)J in the ‘±’ cells via (5.3a), satisfy the heat equation (5.3bl away 
from the cell corners. Using this and the no-flux boundary conditions (3.7) makes it 
possible to write 


$+(C, 2-) = (?{+ +Jf_)$+(C,0+), 
$+(C, 6 -) = (?{++Jf_)S+(C,4+), 


(B.la) 

(B.lb) 

where 34^ are the ‘time-2’ heat-flow maps defined in ( |5.6[ ). Continuity of 4)+ across 
the half-cells implies that 


(B.lc) $+{(, 4-) = 3f_$+(C, 2 +) + 3f+$^(C, 2+), 

(B.ld) $+(C, 0-) = 3f_$+(C, 6 +) + 3{+$L(C, 6 +), 
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where and are the solutions inside the neighbouring ’ cells, located on the 
left and the right of the “+’ cell, respectively. Using the alternating symmetry (3.121 
and the definition (5.3a) further gives 

$«((,2+) =e-"«$+(C,6+), 


(B.le) 

(B.lf) 


$^((,2+) = e’"'^$+(C,2+). 


Employing the jump condition (5.4) to eliminate the upstream functions (at a = 
0“, 2“, 4“ and 6“) in favour of the downstream ones (at cr = 0+, 2+, 4+ and 6+) 
reduces (B.l) to the eigenvalue problem 


(B.2a) 

where 

(B.2b) 


(16Pe)*/2$(C) = (/C$)(C) 


$ = 


/$+(C,o+)\ 

$+(C,2+) 

$+(C,4+) 

V$+(c,6+)y 


K = 


0 

L++L_ 

0 

0 


0 

L+ 

0 


0 

0 

0 

L+ +L_ 


0 

e-’"«L_ 

0 


and L± = 


Appendix C. Matching of Regimes II and III. In this section, we derive 
the form of / in subregime Ilb and show that it matches the corresponding expression 
in subregime Ilia. We seek an approximation to the principal eigenvalue A of Ki - 
or more accurately to log A from which /o is inferred - in the limit g 3> 1 and hence 
/o ^ 1. It turns out to be advantageous to consider 


(C.l) 

noting the approximations 

/$(C,o+)\ 

(C.2) #(C) = 


(/C2$)(C) = A2$(C), 


<1>(C,2+) 

0 

V 0 


+ 0(e-"«) and A2(g,/o) = e"«A(/o) + 0(1). 


Here, A(/o) is the principal eigenvalue of the integral equation 

(C.3) A(/o)<^(0 = (a_(^)(C) + (J+¥^)(C), with = L_L±. 

Note that the above simplification is possible because L_L+ is the adjoint of L+L_ 
(and thus they share the same eigenvalues). 

The asymptotic behaviour of A(/o) for /o 3> 1 is obtained by introducing the 
1 /2 

rescaling ( = f^' z into (C.3). It now becomes natural to employ a WKB expansion 
for the principal eigenfunction so that, at leading order, (p{z) ^ exp(/oS'(z))A(z), 
where S{z) and A(z) remain to be determined. Thus, 


(C.4) (J±^) (z) 


ffo+l 

Jo _ 

Stt 


dzi / dz 2 exp{fo{h±{z,zi,Z 2 ) + S{z 2 ))) A{z 2 ), 
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where, from the definition of in (5.6) 
(C.5) 


h±{z,zi,Z2) = log(2:2:i) - (z =F zi) - {zi + -^2) / 8 . 


The contribution of {0-ip){z) to (C.3) is subdominant with respect to 
because, when /o 3> 1, exp(/oh_) is exponentially smaller than exp(/o/i+) for all 
Zi,Z 2 > 0. Using Laplace’s method in (C.3) we obtain that 

(C.6) a= lim /o"Mog(A(/o)/o"*) =/i+(z, zi, za) + S'(z 2 ) - ^(z), 

/o-^-oc 

where zi(z) and Z 2 (z) are determined by the saddle-point conditions dzih+ = 0 and 
= —5"(z 2). The constant a governs the asymptotics of log A(/o) and hence of 

log A. 

Now, the right-hand side of ( ra can be obtained without knowledge of S{z) if 
it is evaluated at the solution Z of Z 2 {Z) = Z. We now obtain expressions for zi{Z) 
and Z. Note first that differentiation of ( |C.6 ) with respect to z gives dzh+ = S'{z). 
Combining this with the saddle-point conditions leads to 


(C.7) 


dh^ 

dzi 


(Z, zi(Z), Z) = 0 and 


dh^ 

dz 


i9h_| 

dz2 J 


1 


(Z,zi(Z),Z)=0. 


Using the explicit form of in (C.5), these equations are readily solved to find that 
Z = zi(Z) = Z 2 {Z) = v^, whence 


(C.8) 


o = log 2 - 1 and log A(/o) - /o log(2/o/e). 


Employing the latter expression into ( |C.2 ) provides an expression for 2 log A that we 
use inside (5.7) to obtain that /o(log(16Pe) — a) = 7rq -|- /olog/o -I- 0(1). It is now 


relatively straightforward to deduce that 

—nq 


(C.9) 


/o 


W,n (-^g(8ePe)-i) 


for 1 <C g <C Pe, 


where Wm denotes the second real branch of the Lambert W function (see e.g., [T]). 
The upper bound in (C.9) corresponds to the upper value of q for which /o remains 
a non-decreasing function of q. 
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